Loading Libraries¶
In [2]:
library("tidyverse")
library("GGally")
library("corrplot")
library("reshape2")
library("plotly")
library(MASS)
library(ggplot2) # plot
library(calibrate) # textxy
library(car) #vif, durbin-watson test
library(orcutt) # To fix atuto-correlation issue
library(leaps) #regsubset
library(qpcR) #To calculate PRESS Residual and PRESS Statistic
library(factoextra)# PCA graphs
library(nortest) # for normality testing - Anderson-Darling etc
Increasing the resolution of the plots:¶
In [2]:
options(repr.plot.height=8, repr.plot.res=300)
Loading the Dataset¶
In [3]:
boston_dataset <- read_csv("../datasets/BostonHousing.csv")
Data Exploration And Cleaning¶
Glance of the Dataset¶
In [4]:
boston_dataset <- boston_dataset %>%
dplyr::select(MEDV, everything())
head(boston_dataset)
Features of the data:¶
In [5]:
names(boston_dataset)
Following information we got from the official site: https://archive.ics.uci.edu/ml/machine-learning-databases/housing
There are 14 attributes in each case of the dataset. They are:
- CRIM - per capita crime rate by town
- ZN - proportion of residential land zoned for lots over 25,000 sq.ft.
- INDUS - proportion of non-retail business acres per town.
- CHAS - Charles River dummy variable (1 if tract bounds river; 0 otherwise)
- NOX - nitric oxides concentration (parts per 10 million)
- RM - average number of rooms per dwelling
- AGE - proportion of owner-occupied units built prior to 1940
- DIS - weighted distances to five Boston employment centres
- RAD - index of accessibility to radial highways
- TAX - full-value property-tax rate per $\$10,000$
- PTRATIO - pupil-teacher ratio by town
- B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
- LSTAT - %lower status of the population
- MEDV - Median value of owner-occupied homes in $1000's
MEDV - Price of the houses is our target variable.
Dimension of the data:¶
In [6]:
dim(boston_dataset)
So our data has 506 observations of 14 columns.
Datatypes of columns:¶
In [7]:
str(boston_dataset)
We can see that all the features are already in numerical double formats.
Summary of the data:¶
In [8]:
summary(boston_dataset)
In the above summary we can see the helpful statistics such as min, max, mean, median and quartiles for all of the columns.
Data Cleaning¶
Checking Missing Values:¶
In [7]:
colSums(is.na(boston_dataset))
Inference: There is no NULL or missing values in our data. We will further investigate upon this.
Data Visualization¶
Pairwise features relation¶
In [25]:
options(repr.plot.height=12)
ggpairs(boston_dataset)
In the above ggpairs plot we can see the correlation, denisty curve and the scatterplot between pairs simultaenously. Below we will see more plot specifically focused on the correlation and density curves then derive inferences.
Inference:
- Most of the variables are not normal.
In [24]:
options(repr.plot.height=8)
corrplot.mixed(cor(boston_dataset), lower="number", upper="circle",tl.pos = "lt")
Inference: From the above data we can see the correlation between features. We can infer:
- (INDUS, NOX), (AGE, NOX), (TAX, INDUS), (TAX, RAD) and (MEDV, RM) paris are highly positively correlated.
- (DIS, INDUS), (DIS, NOX), (DIS, AGE) and (MEDV, LSTAT) pairs are highly negatively correlated.
Because MEDV is our target variable and we see two highly correlated features with it, we are interested in the jointdensity plot for them. They are given below:
In [12]:
options(repr.plot.height=5, repr.plot.res=200)
# devtools::install_github("WinVector/WVPlots")
suppressMessages(library(WVPlots))
suppressWarnings( ScatterHist( boston_dataset, "MEDV", "LSTAT", title="Joint Density plot of MEDV - LSTAT"))
suppressWarnings( ScatterHist( boston_dataset, "MEDV", "RM", title="Joint Density plot of MEDV - RM"))
Inference:
- MEDV is highly positively correlated with RM
- MEDV is highly negatively correlated with LSTAT.
3D Scatterplot for LSTAT, RM and MEDV:¶
In [8]:
library(rgl)
# Plot
plot3d(
x=boston_dataset$LSTAT, y=boston_dataset$RM, z=boston_dataset$MEDV,
# col = data$CHAS,
type = 's',
radius = .1,
xlab="LSTAT", ylab="RM", zlab="MEDV")
# To display in an R Markdown document:
rglwidget()
# To save to a file:
htmlwidgets::saveWidget(rglwidget(width = 800, height = 800),
file = "3D-Scatterplot-MEDV-LSTAT-RM.html",
libdir = "libs",
selfcontained = FALSE
)
Inference:
- The points in the 3D scatterplot isnt lying on a plane.
- There isnt a proper linear relationship among LSTAT, RM and MEDV.
Checking the outliers and density curve:¶
In [14]:
#setting the resolution of the plots
par(mfrow=c(1,3))
options(repr.plot.height=3,repr.plot.width=15, repr.plot.res=150)
# a function to plot density curve, histogram and boxplot of a feature simultaenously.
powerPlot <- function(x, name)
{
dens <- density(x)
# Histogram
hist(x, probability = TRUE, xlab = name, ylab="", ylim= c(0, max(dens$y)), col = "grey",
axes = FALSE, main = "", cex.lab=2)
# Axis
axis(1, cex.axis=2)
# Density
lines(dens, col = "red", lwd = 1.5)
# Add boxplot
par(new = TRUE)
boxplot(x, horizontal = TRUE, axes = FALSE,
lwd = 2, col = rgb(0, 1, 1, alpha = 0.15), outlier.colour="red")
}
# a loop to run powerPlot() for all of the features
for (x in names(boston_dataset)){
powerPlot( boston_dataset[[x]], x)
}
Inference:
- None of the feature in perfect Normal distribution.
- We can see that CRIM, ZN, RM, B and MEDV have of outliers.
- Also we note that there isnt any specific one transformation that fits all features.
Relation of MEDV against other features:¶
In [15]:
options(repr.plot.height=12, repr.plot.res=200)
rel <- boston_dataset %>%
melt(id.vars = "MEDV") %>%
ggplot(aes(x = value, y = MEDV, colour = variable), na.rm=TRUE) +
geom_point(alpha = 0.7, na.rm=TRUE) +
stat_smooth(aes(colour = "black"), na.rm=TRUE) +
facet_wrap(~variable, scales = "free", ncol = 2) +
labs(x = "Variable Value", y = "Median House Price ($1000s)") +
theme_minimal( base_size=18) +
theme( legend.text = element_text( size=14))
suppressMessages(suppressWarnings( print(rel)))
Inference:
- There is no perfect linear relationship between MEDV and any other feature.
- Most of the relationship is non-linear.
- CHAS is a categorical feature, more concentrated on CHAS=0.
Density Curve of Target Variable (MEDV - Price of the houses)¶
In [32]:
options(repr.plot.height=5, repr.plot.res=200)
boston_dataset %>%
ggplot(aes(MEDV)) +
stat_density(fill="white", col="blue") +
theme_minimal( base_size=18) +
ggtitle("MEDV-Price of houses Density curve") + theme( plot.title = element_text( hjust=0.5))
Inference:
- The target variable is not normal.
Checking Outliers¶
In [33]:
count=0
outliers=c()
for(i in 1:nrow(boston_dataset)){
x <- boston_dataset[i,]
for(col in names(boston_dataset)){
mu=mean(boston_dataset[[col]])
sigma=sd(boston_dataset[[col]])
if(x[[col]]<mu-3*sigma || x[[col]]>mu+3*sigma)
{
count=count+1
}
}
if(count>0){
outliers <- c(outliers, i)
}
count=0
}
cat("No. of Outliers:",length(outliers))
MODELLING ¶
Feature Engineering¶
Backward Elimination¶
In [34]:
## backward elimination
BER=regsubsets(MEDV~.,data=boston_dataset,method="backward")
Modelsummary=cbind(summary(BER)$which,R2=summary(BER)$rsq,SSres=summary(BER)$rss,AdjR2=summary(BER)$adjr2,Cp=summary(BER)$cp,BIC=summary(BER)$bic)
Modelsummary
Significant Features : CRIM , NOX , RM , DIS , RAD , PTRATIO , B , LSTAT
Forward Substitution¶
In [35]:
## forward substituition
FSR=regsubsets(MEDV~.,data=boston_dataset,method="forward")
# summary(FSR)
Modelsummary=cbind(summary(FSR)$which,R2=summary(FSR)$rsq,SSres=summary(FSR)$rss,AdjR2=summary(FSR)$adjr2,Cp=summary(FSR)$cp,BIC=summary(FSR)$bic)
Modelsummary
Significant Features : ZN , CHAS , NOX , RM , DIS , PTRATIO , B , LSTAT
Stepwise Regression¶
In [36]:
#Stepwise Regression
SWR=regsubsets(MEDV~.,data=boston_dataset,method="seqrep")
Modelsummary=cbind(summary(SWR)$which,R2=summary(SWR)$rsq,SSres=summary(SWR)$rss,AdjR2=summary(SWR)$adjr2,Cp=summary(SWR)$cp,BIC=summary(SWR)$bic)
Modelsummary
Significant Features : ZN , CHAS , NOX , RM , DIS , PTRATIO , B , LSTAT
Model-1:¶
We will create a model taking the intersection of features suggested by the above three methods.
The glance of data of the suggested features is given below.
In [5]:
boston_dataset <- boston_dataset %>% dplyr::select(c(MEDV,ZN,CHAS,NOX,RM,DIS,PTRATIO,B,LSTAT))
head(boston_dataset)
Modelling:¶
In [7]:
modelm <- lm(boston_dataset$MEDV ~ ., data =boston_dataset )
summary(modelm)
# RMSE for checking the accuracy of model
cat("RMSE:", sqrt(sum(modelm$residuals^2)/modelm$df))
Model-1
RMSE = 4.84
Adjusted R-Squared = 0.722
Adequecy Checking of Model-1¶
Checking the multi-collinearity in model-1:¶
In [11]:
vif(modelm)
Inference:
- Since all values are well under 5, thus we can conclude that the features are not multicollinear.
Checking Influence Points (Model-1)¶
In [8]:
## Checking for Influence Points
inflm.boston <- influence.measures(modelm)
inf_pts<-which(apply(inflm.boston$is.inf, 1, any))
# which observations 'are' influential
dim(summary(inflm.boston)) # only these
Inference:
- There are 54 potentially influental points in our data.
Since this can affect our model adversely we check which influence points lies outside 2-sigma limits and remove those points.
Influence points that are outliers:¶
In [41]:
#boston_dataset<-boston_dataset[-inf_pts,]
count=0
outliers=c()
for(i in inf_pts){
x <- boston_dataset[i,]
for(col in names(boston_dataset)){
mu=mean(boston_dataset[[col]])
sigma=sd(boston_dataset[[col]])
if(x[[col]]<mu-2*sigma || x[[col]]>mu+2*sigma)
{
count=count+1
}
}
if(count>0){
outliers <- c(outliers, i)
}
count=0
}
length(outliers)
Inference:
- 51 out of 54 influence points lie outside 2-sigma limits.
Removing all 54 points resulted in better model.
In [9]:
# removing outlier-influence-points
boston_dataset<-boston_dataset[-inf_pts,]
dim(boston_dataset)
Modelling after removing the influental points - Model-1.2:¶
In [10]:
names(boston_dataset)
modelm <- lm(boston_dataset$MEDV ~ ., data =boston_dataset )
summary(modelm)
cat("RMSE:", sqrt(sum(modelm$residuals^2)/modelm$df))
Model-1.2
RMSE = 3.59
Adjusted R-Squared = 0.81
Thus we got considerable improvement in RMSE from 4.84 to 3.59.
Multicollinearity Check¶
In [15]:
vif(modelm)
Inference: There is no multicollinearity.
Normality Check¶
In [21]:
par(mfrow=c(1,3))
options(repr.plot.height=3,repr.plot.width=15, repr.plot.res=150)
powerPlot <- function(x, name) {
dens <- density(x)
# Histogram
hist(x, probability = TRUE, xlab = name, ylab="", ylim= c(0, max(dens$y)), col = "grey",
axes = FALSE, main = "", cex.lab=2)
# Axis
axis(1, cex.axis=2)
# Density
lines(dens, col = "red", lwd = 1.5)
# Add boxplot
par(new = TRUE)
boxplot(x, horizontal = TRUE, axes = FALSE,
lwd = 2, col = rgb(0, 1, 1, alpha = 0.15), outlier.colour="red")
}
# a loop to run powerPlot() for all of the features
for (x in names(boston_dataset)){
powerPlot( boston_dataset[[x]], x)
}
Inference:
- The data are not perfectly normal.
Checking Normality of Residuals (Model-1.2):¶
- QQ Plot
- Shapiro Test
- Anderson-Darling Test
In [103]:
??qqnorm
In [22]:
# Test-1. QQ Plot
options(repr.plot.height=8, repr.plot.res=300)
ors <- rstandard(modelm)
qqnorm(ors, pch = 1, frame = FALSE, cex.lab=1.5, cex.axis=1.5, main="Normal QQ Plot of Residuals", cex.main=2)
qqline(ors, col = "steelblue", lwd = 2)
In [50]:
# Test-2
shapiro.test(rstandard(modelm))
In [51]:
# Test-3: Anerson Darling Test
ad.test(rstandard(modelm))
Inference: From the above tests we can infer following:
- The residuals are not normal.
We estimate this result due to the fact that the target variable in the data wasnt distributed normally.
Linear Relation between Residuals and Xs¶
In [52]:
## linear relation between residuals and X's
par(mfrow=c(3,3))
for (x in names(boston_dataset)){
plot(boston_dataset[[x]],rstudent(modelm),xlab = x,ylab = "rst(ti)", cex.lab=2, cex.axis=2)}
Inference:
- The X values are not linearly related to Y.
Principal Component Regression Analysis¶
In [23]:
# GETTING Principal Components
boston_x <- boston_dataset[,c(-1)]
#Scale the data
res.pca <- prcomp(boston_x, retx=TRUE,center = FALSE, scale = TRUE)
Summary of PCAs:¶
In [24]:
summary(res.pca)
Eigen values or equivalently variance of Principal Components: (The same information is available in the above summary)
In [25]:
eigval<- get_eigenvalue(res.pca)
eigval
Scree Plot¶
In [26]:
options( repr.plot.height=5)
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0,100), ggtheme = theme_minimal(base_size=18)+ theme( plot.title = element_text( hjust=0.5)))
Inference:
- PC-1 explains the most proportion of variance.
- The elbow occurs at PC-2.
Below we will look further into PC-1 and PC-2.
Checking contributions of variables to Principal Components¶
In [27]:
library(factoextra)
options(repr.plot.height=5, repr.plot.res=300)
# Contributions of variables to PCs
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10, ggtheme = theme_minimal(base_size=18)+ theme( plot.title = element_text( hjust=0.5)))
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10, ggtheme = theme_minimal(base_size=18)+ theme( plot.title = element_text( hjust=0.5)))
Inference:
- PC-1 is contributed largely by multiple variables.
- PC-2 is largely contributed mostly by CHAS.
The below chart shows the same thing.
In [140]:
fviz_pca_var(res.pca, col.var = "blue", ggtheme = theme_minimal(base_size=18))
Model with PCA1 suggested features data: Model-2¶
In [19]:
##model with PCA1 suggested features data.
model_pca2 = lm(MEDV~RM+PTRATIO+DIS+LSTAT,data=boston_dataset)
summary(model_pca2)
cat("Adjusted R squared:",summary(model_pca2)$adj.r.squared,"\n")
cat("RMSE:", sqrt(sum(model_pca2$residuals^2)/model_pca2$df))
Model-2
RMSE = 3.84
Adjusted R-squared = 0.79
lm() Model¶
The lm() model is built based upon the suggested features of all of the previous methods. The influental outliers also have been removed in previous section.
Modelling:¶
In [18]:
lm_model = lm(MEDV~., data=boston_dataset)
summary(lm_model)
cat("Adjusted R squared:",summary(lm_model)$adj.r.squared,"\n")
cat("RMSE:", sqrt(sum(lm_model$residuals^2)/lm_model$df))
FInal Model
RMSE = 3.59
Adjusted R-Squared = 0.81
Above is the best value we got for RMSE and Adjusted R-squared among all of our models.
Predicting using the model:¶
In [28]:
options( repr.plot.height=8)
pred<- predict(lm_model, boston_dataset)
act_vs_pred<- data.frame(actual= boston_dataset$MEDV, predicted= pred)
ggplot(act_vs_pred, aes(x= predicted, y=actual))+
geom_point()+
geom_abline(intercept=0, slope=1)+
theme_minimal(base_size=18) + theme( plot.title = element_text( hjust=0.5)) +
labs(x= 'Predicted Values', y= "Actual Values", title= "Predicted vs Actual Values")
Inference:
- There is the linearity between predicted values and the observed values.
Checking the normality of residuals:¶
In [29]:
options(repr.plot.res=300, repr.plot.height=8)
res <- as.data.frame(residuals(lm_model))
ggplot(res,aes(residuals(lm_model)))+
geom_histogram(fill='brown',alpha=0.5, binwidth=1) +
theme_minimal( base_size=18) + theme( plot.title = element_text( hjust=0.5)) + labs(title="Distribution of Residuals")
Inference:
- The residuals does not follow the exact normal distribution.
- But we take it to be normal within accepted margins.
In [30]:
options(repr.plot.height=12,repr.plot.res=300)
par(mfrow = c(2,2))
plot(lm_model, cex.lab=1.5, cex.axis=1.5)
CV Model¶
Partitioning training and test dataset:¶
In [11]:
suppressMessages(library(caret))
set.seed(998)
trainPartitionRows <- createDataPartition(boston_dataset$MEDV, p = .90, list = FALSE)
cat("Total no of observations = ", nrow(boston_dataset),"\n")
cat("# for training = ", nrow(trainPartitionRows),"\n")
cat("# for testing = ", nrow(boston_dataset) - nrow(trainPartitionRows),"\n")
In [12]:
trainDataset <- boston_dataset[ trainPartitionRows,]
testDataset <- boston_dataset[-trainPartitionRows,]
cat("\nHead of TRAIN DATASET:")
head(trainDataset)
cat("\nHead of TEST DATASET:")
head(testDataset)
Performing Cross-Validation:¶
In [13]:
fitControl <- trainControl(## 10-fold CV
method = "cv", # repeatedcv
number = 10
## repeated ten times
#repeats = 10,
)
CVModel <- train(MEDV ~ .,
data = trainDataset,
method = "lm",
trControl = fitControl,
## This last option is actually one
## for gbm() that passes through
#verbose = TRUE
)
CVModel
summary(CVModel)
CV Model
RMSE = 3.55
Adjusted R-squared = 0.82
Mean Absolute Error (MAE) = 2.69
Validation on training data for CV-Model:¶
In [31]:
#validation on training data
pred<- predict(CVModel, trainDataset)
validation<- data.frame(actual= trainDataset$MEDV, predicted= pred)
ggplot(validation, aes(x= predicted, y=actual))+
geom_point()+
geom_abline(intercept=0, slope=1)+
theme_minimal(base_size=18) + theme( plot.title = element_text( hjust=0.5)) +
labs(x= 'Predicted Values', y= "Actual Values", title= "Predicted vs Actual Values")
Inference:
- Their is the linear relation between observed and predicted values.
Summary of Prediction:¶
In [29]:
residuals<- (validation$actual - validation$predicted)^2
RMSE <- sqrt(mean(residuals))
summary(pred)
cat("RMSE on train data:",RMSE)
Validation on Test data:¶
In [35]:
#TESTING
pred1<- predict(CVModel, newdata= testDataset)
test<- data.frame(actual= testDataset$MEDV, predicted= pred1)
ggplot(test, aes(x= predicted, y=actual))+
geom_point()+
geom_abline(intercept=0, slope=1)+
theme_minimal(base_size=18) + theme( plot.title = element_text( hjust=0.5)) +
labs(x= 'Predicted Values', y= "Actual Values", title= "Predicted vs Actual Values")
Summary of validation on test data:¶
In [32]:
summary(pred1)
residuals1<- (test$actual - test$predicted)^2
RMSE1 <- sqrt(mean(residuals1))
cat("RMSE on test data:",RMSE1)
Inference:
RMSE on test data: 3.94630
RMSE on train data: 3.513884
- The difference in RMSE between train and test data is ~0.4
Training perfomance vs Testing performance¶
In [163]:
n=nrow(trainDataset)
# create empty data frame
learnCurve <- data.frame(m=double(n),
trainRMSE=double(n),
cvRMSE = double(n))
# test data response feature
testY <- testDataset$MEDV
# Run algorithms using 10-fold cross validation with 3 repeats
trainControl <- trainControl(method="repeatedcv", number=10, repeats=3)
metric <- "RMSE"
# loop over training examples
for (i in 137:nrow(trainDataset)) {
learnCurve$m[i] <- i
# train learning algorithm with size i
fit.lm <- train(MEDV~.,
data=trainDataset[1:i,],
method="lm",
metric="RMSE",
#preProc=c("center", "scale"),
trControl=trainControl)
learnCurve$trainRMSE[i] <- fit.lm$results$RMSE
# use trained parameters to predict on test data
prediction <- predict(fit.lm, newdata = testDataset[,-1])
rmse <- postResample(prediction, testY)
learnCurve$cvRMSE[i] <- rmse[1]
}
In [164]:
plot(log(learnCurve$trainRMSE),type = "o",col = "red",xlab = "Training set size",
ylab = "Error (RMSE)", main = "Linear Model Learning Curve", ylim=c(1,1.50), xlim=c(136,410), cex.axis=1.5, cex.lab=1.5)
lines(log(learnCurve$cvRMSE), type = "o", col = "blue")
legend('topright', c("Train error", "Test error"), lty = c(1,1), lwd = c(2.5, 2.5),
col = c("red", "blue"))
Saving Model¶
Since the best performing model is CV model thus we take same as the CV model.
In [13]:
final_model <- CVModel
saveRDS(final_model, "./model.rds")
In [204]:
summary(final_model)
Final Model(CV Model)
RMSE = 3.55
Adjusted R-squared = 0.82
Mean Absolute Error (MAE) = 2.69